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Abstract 



> 

We present a theoretical background for the data analysis of the gravitational-wave signals from 

■ spinning neutron stars for Earth-based laser interferometric detectors. We introduce a detailed model 
' of the signal including both the frequency and the amplitude modulations. We include the effects of 
, the intrinsic frequency changes and the modulation of the frequency at the detector due to the Earth 

■ motion. We estimate the effects of the star's proper motion and of relativistic corrections. Moreover 
we consider a signal consisting of two components corresponding to a frequency / and twice that 

CJ> , frequency. From the maximum likelihood principle we derive the detection statistics for the signal 

and we calculate the probability density function of the statistics. We obtain the data analysis 
I I procedure to detect the signal and to estimate its parameters. We show that for optimal detection 

bJQ. of the amplitude modulated signal we need four linear filters instead of one linear filter needed for a 

constant amplitude signal. Searching for the doubled frequency signal increases further the number of 
linear filters by a factor of two. We indicate how the fast Fourier transform algorithm and resampling 

■ methods commonly proposed in the analysis of periodic signals can be used to calculate the detection 
^ ' statistics for our signal. We find that the probability density function of the detection statistics 

is determined by one parameter: the optimal signal-to-noise ratio. We study the signal-to-noise 
ratio by means of the Monte Carlo method for all long-arm interferometers that are currently under 
construction. We show how our analysis can be extended to perform a joint search for periodic signals 
by a network of detectors and we perform Monte Carlo study of the signal-to-noise ratio for a network 
of detectors. 



1 Introduction 

Spinning neutron stars are one of the primary candidate sources of gravitational waves for long-arm laser 
interferometric detectors (|^, see for a review). Detectors with a sufRcent sensitivity to see strong 
neutron star sources anywhere in the Galaxy will be taking data within two or three years ^, |[ ^. 
A rotating body, perfectly symmetric about its rotation axis does not emit gravitational waves. If the 
spinning neutron star is to emit gravitational waves over extended periods of time, it must have some kind 
of long-lived asymmetry. Several mechanisms have been given for such an asymmetry to arise ^, ^ . 
During the crystalization period the crust of the neutron star may develop deviations from axisymmetry 
that will be supported by anisotropic stresses in the solid crust The strong magnetic field present 
in the neutron star may not be aligned with the rotation axis and consequently the distortion produced 
by the magnetic pressure results in the neutron star being asymmetric Also the rotation axis may 
not coincide with a principal axis of the star's moment of inertia tensor. Then the star will precess and 
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emit gravitational waves §, There are other mechanisms that can produce gravitational waves from 
neutron stars. Accretion of matter on a neutron star can drive it into a nonaxisymmetric configuration 
and power steady radiation with a considerable amplitude. This mechanism has been pointed out by 
Wagoner ([0, see also [|l2|). It applies to a certain class of neutron stars, including accreting stars in 
binary systems that have been spun up to the first instability point of the Chandrasekhar-Friedman- 
Schutz (CFS) instability [|l3[ Recently Andersson suggested a similar instability in r— modes 
of rotating relativistic stars. The effectiveness of these instabilities depends on the viscocity of the star 
which in turn is determined by the temperature of the star . 

This paper initiates a series of papers where theoretical problems of data analysis of gravitational- 
wave signals from spinning neutron stars are considered, independently of the mechanisms generating the 
waves. 

The data analysis of monochromatic signals for intcrferomctric antennae was investigated by one of us 
p5| . A search strategy for such signals was proposed and the computing power required estimated. The 
basic method to detect periodic signals is to Fourier analyse the data, and an efficient computational tool 
is the fast Fourier transform. The main problem is that to do the search one has to take into account the 
modulation of the signal due to the Earth's motion relative to the solar system barycenter. If the position 
of the source on the sky is unknown this introduces two additional parameters in the signal and this vastly 
increases the computational time to do the search. It is clear that the main limit on the sensitivity of 
such a search will be the available computing power. Variants of the proposed search strategy have been 
implemented with test data from the prototype detectors where the search was carried out only over a 
limited region of the parameter space . 

The problem of computational requirements has recently been reconsidered by Brady et al. 
They realized that in the model of the signal the effect of the intrinsic frequency modulation due to 
spin-down or spin-up of the neutron star needs to be considered. This increases the parameter space and 
consequently the computational power required to search all the parameter space. Assuming access to 
teraflops computing power it was shown that coherent integration times will be limited to days for an 
all-sky search for young, rapidly spinning stars and to weeks for more directed searches. A simplified 
model of the signal where modulation due to diurnal rotation of the Earth was neglected has also been 
examined by one of us [ p2[ and the computational requirements to do the search were estimated. 

In this series of papers we consider a more general model of the signal than in the work cited above. We 
take into account not only the modulation of the phase of the signal but also the amplitude modulation. 
Moreover we consider a signal consisting of two components corresponding to a frequency / and twice 
that frequency. In general neither of the components is dominant. 

In this work, which is Paper I of the series, we introduce the signal and we derive an optimal data 
analysis procedure for its detection. In Paper II we examine the accuracy of estimation of the parameters 
of the signal that can be achieved with the optimal analysis. In Paper III we examine in detail the 
characteristics of the detection statistics derived in Paper I and the computational power required to 
calculate it. In Paper IV we investigate the least-squares method to estimate astrophysically interesting 
parameters of the signal from the estimators of the amplitudes derived in Paper I. 

The plan of this paper is as follows. In Section 2 we derive a general formula for the response of a laser 
interferometer to our two component signal including both the phase and the amplitude modulation. In 
Section 3 from the maximum likelihood principle we derive the data analysis procedure to detect the signal 
introduced in Section 2 and to estimate its parameters. We obtain the basic probability density functions 
of the detection statistics. We show that probability of detection is determined by one parameter: the 
optimal signal-to-noise ratio. We study this quantity by means of the Monte Carlo simulations for all the 
interferometric detectors that are currently under construction. We conclude Section 3 by showing how 
one can take advantage of the speed of the FFT algorithm to evaluate efficiently our detection statistics. 
This involves application of the resampling techniques proposed earlier for the case of a simpler signal 
model |2l[. In Section 4 we show that our analysis can easily be extended to networks of detectors 
and wc perform a Monte Carlo study of the signal-to- noise ratio for the networks. In Appendix A we 
discuss the model of the phase of the gravitational-wave signal and in particular we estimate the effect of 
the proper motion of the neutron star and of relativistic corrections. In Appendix B we give the general 
analytic formula for the optimal signal-to-noise ratio. In Appendix C we present an additional study by 
means of the Monte Carlo simulations of the signal-to-noise ratios for the individual components of the 
spinning neutron star signal. 
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2 Noise-free response of the interferometric detector 



2.1 Beam-pattern functions 

The response of a laser interferometric detector to a weak plane gravitational wave in the long wavelength 
approximation (i.e. when the size of the detector is much smaller than the reduced wavelength A/(27r) of 
the wave) is well known (see, e.g., and Section IIA of |Q and references therein). The dimensionless 
detector's response function h is defined as the difference between the wave induced relative length changes 
of the two interferometer arms and can be computed from the formula (cf. Eq. (5) of [p3) 
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/i(i) = -ni- H{t)ni ~-n2- H{t)n2 
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(1) 



where rii and n2 denote the unit vectors parallel to the arm number 1 and 2, respectively (the order of 
arms if defined such that the vector ni x n2 points outwards from the surface of the Earth), H is the 
3-dimensional matrix of the spatial metric perturbation produced by the wave in the proper reference 
frame of the detector, and a dot stands for the standard scalar product in the 3-dimensional Cartesian 
space. The matrix H is given by ^ 

H{t) = M{t)H{t)M{tf , (2) 

where M is the 3-dimensional orthogonal matrix of transformation from the wave Cartesian coordinates 
{xwj Uw, Zw) to the Cartesian coordinates {xd, yd, -^d) in the detector's proper reference frame (the defini- 
tion of these coordinates is given below), T denotes matrix transposition. In the wave coordinate system 
the gravitational wave travels in the +Zw direction. In this frame the matrix H has the form 

/ h+{t) h^i) \ 
H{t)= h^t) -h+{t) , (3) 
\ / 

where the functions h+ and describe two independent wave's polarizations. Collecting Eqs. (|l|)-(||) 
together one can see that the response function /i is a linear combination of the functions and h-^: 

h{t)^F+{t)h+{t)+FAt)hy,{t), (4) 

where F+ and F-^ are called the beam-pattern functions. 

Because of the diurnal motion of the Earth the beam-patterns and Fx are periodic functions of 
time with a period equal to one sidereal day. We want now to extract explicitly this time dependence as 
well as to express F^ and -Fx as functions of right ascension a and declination 6 of the gravitational- wave 
source and polarization angle 4' (the angles a, S, and V' determine the orientation of the wave reference 
frame with respect to the celestial sphere reference frame defined below) . Our treatment partially follows 
that of Section 5 of ||^ . We represent the matrix M of Eq. (H) as 

M = AhAhMl , (5) 

where Mi is the matrix of transformation from wave to celestial sphere frame coordinates, M2 is the 
matrix of transformation from celestial coordinates to cardinal coordinates and is the matrix of 
transformation from cardinal coordinates to detector proper reference frame coordinates. In celestial 
sphere coordinates the z axis coincides with the Earth's rotation axis and points toward the North pole, 
the X and y axes lie in the Earth's equatorial plane, and the x axis points toward the vernal point. In 
cardinal coordinates the (x, y) plane is tangent to the surface of the Earth at detector's location with 
X axis in the North-South direction and y axis in the West-East direction, the z cardinal axis is along 
the Earth's radius pointing toward zenith. In detector coordinates the z axis coincides with the z axis 
of cardinal coordinates and the x axis is along the first interferometer arm (then the y axis is along the 
second arm if the arms are at a right angle) . Under the above conventions the matrices Mi , M2 , and Afa 
are as follows (matrices Mi and M2 given below coincide with matrices A and B from Ref. cf. Eqs. 
(52) and (60) of [|) 

sin a cos ■0 — cos a sin 5 sin ■0 — cos a cos ip — sin a sin S sin ip cos 6 sin ip 
Ml = I — sin a sin — cos a sin (5 cos ^0 cos a sin -0 — sin a sin 5 cos "0 cos (5 cos 0" | , (6) 
— cos a cos S — sin a cos S — sin 6 
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^Irt) sin X Sm[(j)r + flrt) 
M2 = I - Sin{(j)r + rirt) COs{(j)r + flrt) | , (7) 

ilrt) COS A Sm{(j)r + i^rt) 

- sin (7 + C/2) cos (7 + C/2) \ 
M3 = I -cos (7 + C/2) -sin (7 + C/2) . (8) 
1/ 

In Eq. A is the latitude of the detector's site, Sir is the rotational angular velocity of the Earth, and 
is a deterministic phase which defines the position of the Earth in its diurnal motion at t = (the 
sum (j)r + ^rt coincides with the local sidereal time of the detector's site, i.e. with the angle between the 
local meridian and the vernal point). In Eq. (^) 7 determines the orientation of the detector's arms with 
respect to local geographical directions: 7 is measured counter-clockwise from East to the bisector of the 
interferometer arms, and C is the angle between the interferometer arms. The vectors rii and 112 from 
Eq. (|l|) in the detector's reference frame have coordinates 

ni = (1,0,0), n2 = (cosC,sinC,0). (9) 

The values of the angles A, 7, and the longitudes L (measured positively westwards) for different 
detectors can be found in Table 1 [E5|. 



detector 


A (degrees) 


L (degrees) 


7 (degrees) 


C (degrees) 


GEO600 


52.25 


-9.81 


68.775 


94.33 


LIGO Hanford 


46.45 


119.41 


171.8 


90 


LIGO Livingston 


30.56 


90.77 


243.0 


90 


VIRGO 


43.63 


-10.5 


116.5 


90 


TAMA300 


35.68 


-139.54 


225.0 


90 



Table 1: Positions and orientations of detectors. 



To find the explicit formula for and we have to combine Eqs. (|i|) 
manipulations we arrive at the expressions: 

F+{t) = si-aC [a(t) cos 2tp + b{t)sin2ip] , 
Fx(t) = sinC [b{t) cos 2ip - a(t) sin 2'ip], 



After extensive algebraic 

(10) 
(11) 



where 



ait) 



16 
1 

~4 
1 

+ 4 
1 

~2 
3 



sin27(3 — cos2A)(3 — cos 2(5) cos[2(q! — 0^ — flrt)] 
cos27sinA(3 — cos 2(5) sin[2(Q; ~ (j)r — Sl,-t)] 



sin 27 sin 2A sin 25 cos[a — (pr — ^rt] 
cos 27 cos A sin 26 sin [a — 4>r — 



sin 27 cos^ A cos^ 5, 



cos 27 sin A sin (5 cos [2 ((3! — (pr — ^rt)] 

+ - sin 27(3 — cos 2A) sin 5 sin[2((3: — (^,. — Q.rt)\ 

+ cos 27 cos A cos 5 cos [a — cpr — ^rt] 

+ — sin 27 sin 2A cos 5 sin[a — (p^ — ilrt]. 



(12) 



(13) 



By means of Eqs. (^^-([1^) the beam-pattern functions can be computed directly for any instant of time. 

Equivalent explicit formulae for the beam-pattern functions F+ and Fx (for the case C = '^/'^) can be 
found in Ref. [ p6) where different angles describing the position of the gravitational-wave source in the 
sky and the orientation of the detector on the Earth are used. Also for the case ( = 7t/2 the functions a 
and b can be found in Ref. |l^ , where still another set of angles is used |^ . 
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2.2 The phase of the gravitational-wave signal 

In Appendix A we derive the time dependence of the phase of the gravitational-wave signal observed at 
the detector's location. We consider the significance of the corrections due to the motion of both the 
detector and the neutron star with respect to the the solar system barycenter (SSB) reference frame as 
well as the importance of relativistic corrections. On the basis of the discussion presented in Appendix 
A we adopt the following model of the phase of the gravitational-wave signal: 

^' (fc) fk+l ^ ^ (fc) fk 

nt) = $0 + 2. E foj—^ + -no • r,it) ^ /„-, (14) 

fc=0 ^ k=0 

where /q is the fcth time derivative of the instantaneous frequency evaluated at t = at the SSB, no is 
the constant unit vector in the direction of the star in the SSB reference frame and r^i is the position 
vector of the detector in that frame. 

The signal analysis presented in the remaining part of the paper does not depend on the number s of 
the spindown parameters and therefore we keep s unspecified. 

We associate a coordinate system with the SSB reference frame. The x axis of the system is parallel 
to the X axis of the celestial sphere coordinate system, the z axis is perpendicular to the ecliptic and 
coincides with the orbital angular momentum vector of the Earth. In that system the unit vector no 
pointing towards the star has the components 

\ / cos a cos (5 

cose sine sinacos(5 | , (15) 

— sin e cos e J \ sin 6 

where e is the angle between ecliptic and the Earth's equator. The position vector rd of the detector has 
in this coordinate system the components 

cos A cos (0r + ^rt) 

= Res I sin (00 + ilo^) | + Re | cose sine | | cos A sin (0^ riri) | , (16) 

sin A 

where Res = 1 AU is the mean distance from the Earth's center to the SSB, Re is the mean radius of 
the Earth, flo is the mean orbital angular velocity of the Earth and (po is a deterministic phase which 
defines the position of the Earth in its orbital motion at i = 0. We recall that we neglect the eccentricity 
of the Earth's orbit and the motion of the Earth around the Earth-Moon barycenter. 
Substituting Eqs. (|l|) and (|l|) into Eq. (^ one gets 

*(<) = $o + $(t), (17) 

(fc) ^fc+l 

*W = ^'^g^"(fcTT)! 

H {Res [cos a cos S cos (0o + ^ot) + (cos e sin a cos (5 -f sin e sin S) sin ((j)o + Oot)] 

c 

+ Re [sin X sin S + COS X cos S cos (a — 4>r — ^rt)]} / /oTT- (18) 

fc! 

k=0 

2.3 Wave polarization functions 

We use the following two-component model of the gravitational-wave signal: 

h{t)^hi{t)+h2{t), (19) 




where 



/ii(t) -F+(<)/ii+(i)+Fx(<)/iix(<), h2it) = F+it)h2+{t) + F^it)h2^{t), (20) 
hi+(t) = ■|/ioSin26lsin2tcos*(t), /i2+(t) = jft-o sin^ 0(1 + cos^ t) cos 25'(t), (21) 
hixit) = jho sin 20 sin t sin (t) , ft,2 x (t) = ho sin^ 9 cos t sin 2^1' (t) . (22) 
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The beam-pattern functions F+, Fx are given by Eqs. (|lO|)-([r^) and the phase ^ is given by Eqs. ( |l7| ) 
and 

The model of the signal defined by Eqs. (p^)-(p^) represents the quadrupole gravitational wave that 
is emitted by a freely precessing axisymmetric star. The angle 9, called the wobble angle, is the angle 
between the total angular momentum vector of the star and the star's axis of symmetry and t is the angle 
between the total angular momentum vector of the star and the direction from the star to the Earth. 
The amplitude ho is given by 

ho = — ^, (23) 

where / is the sum of the frequency of rotation of the star and the frequency of precession, / is the 
moment of inertia with respect to the rotation axis, e is the poloidal ellipticity of the star and r is the 
distance to the star. For small wobble angle the signal hi is dominant. Details of the model can be 
found in When = 7r/2 the hi component vanishes. For this special case the /i2 component is the 
quadrupole wave from a triaxial ellipsoid rotating about a principal axis with frequency /. In this case 
the amplitude ho is again given by Eq. (p3|) except that e is now the ellipticity of the star defined by 



(24) 



where Ii and I2 are the moments of inertia of the star with respect to the principal axes orthogonal to 
the rotation axis. This model was considered in |2^. 

Replacing the physical constants in Eq. (E3) by their numerical values results in 



where 



lo = 4.23 X IQ-'^^do ( — - — ^ 
V 100 Hz/ 

~ (iCps) ( 10*5 g err 



(25) 



(26) 



By means of Eqs. (0) and (|ll]) the signal described by Eqs. (|l|)-(||) can be written in the form 



h{t) =^Alihu(t) + ^A2^ h2i{t), 



(27) 



where the eight amplitudes An and A2i are given by 



All 


= ho sin C sin 29 


A12 


— ho sin (^sin29 


Ai3 


= ho sin C sin 29 


An 


— ho sin C sin 29 


A21 


= ho sin C sin^ 9 


A22 


= ho sin C sin^ 9 


A23 


= ho sin ( sin^ 9 


A2i 


— ho sin C sin^ 9 



— sin 2l cos 2-0 cos $0 sin l sin 2^Jj sin <i>o 

8 ^ ^ 4 

— sin L cos 2?/' sin $0 H — sin 2t sin 2ip cos $0 
4 8 

— - sin 2l cos 2^1; sin ^0 — sin t sin 2-0 cos $0 
8 4 

— sin L cos 20 cos $0 ^ g sin 2i sin 2-0 sin <i>o 
-(1 + cos^ i) cos 20 cos 2$o — cos t sin 2-0 sin2<i>o 
i(l + cos^ t) sin20;cos2'i>o + cosiCOs2i/) sin2<i>o 
— -(1 + cos^ l) COS 20 sin 2$o — cos t sin 2-0 cos 2$o 
— -(1 + cos^ l) sin 2-0 sin 2<i>o + cos t cos 20 cos 2$o 



(28) 
(29) 
(30) 
(31) 
(32) 
(33) 
(34) 
(35) 



The amplitudes An and A2i depend on the parameters ho, 9, ip, t, and $0. They also depend on the angle 
Q. The time dependent functions hu have the form 



hii = a{t) cos l<^{t), hi2 = b{t) cosl^{t), 
hi3 = a{t)sml^{t), hi4 ^ 6(i) sin/$(<). 



1 = 1,2, 



(36) 
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where the functions a and b are given by Eqs. and (|T^), respectively, and $ is the phase given by 
Eq. (p^. The modulation amplitudes a and b depend on the right ascension a and the declination S 
of the source (they also depend on the angles A and 7). The phase $ depends on the frequency /o, s 

((c) (fc) 

spindown parameters /g {k = l,...,s), and on the angles a, S. We call parameters /o, /q, a, 5 the 
phase parameters. Moreover the phase $ depends on the latitude A of the detector. The whole signal h 

depends on 8 + s unknown parameters: hg, 0, tp, l, $0i ct, S, /o, /q- 

It is useful to consider the frequency domain characteristics of our gravitational-wave signal. The 
signal consists of two components with carrier frequencies /o and 2/o that are both amplitude and 
phase modulated. The amplitude modulation, determined by functions a and b, splits each of the two 
components into five lines corresponding to frequencies /o — 2/^, fo — fr, fo, fo + fr, fo + '^fr, where 
fr is the frequency of rotation of Earth (/,. ~10~^ Hz) and the same for frequency 2/o. The frequency 
modulation broadens the lines. For the extreme case of the gravitational-wave frequency of 10"^ Hz, the 
spindown age t = 40 years, and the observation time To = 120 days the maximum frequency shifts due 
to the neutron star spindown. Earth orbital motion and Earth diurnal motion are respectively ^--^8 Hz, 
^0.1 Hz, and ~10~'^ Hz. As an example in Figure 1 we have plotted the power spectrum of the noise-free 
response of a detector located near Hannover to the gravitational wave from the Crab pulsar. We took 
only the component /i2 with twice the rotational frequency. We have generated a 24-day long signal. 



3 Optimal filtering for the amplitude modulated signal 
3.1 Maximum likelihood detection 

The signal given by Eq. ( p7| ) will be buried in the noise of a detector. Thus we are faced with the problem 
of detecting the signal and estimating its parameters. A standard method is the method of maximum 
likelihood detection which consists of maximizing the likelihood function A with respect to the parameters 
of the signal. If the maximum of A exceeds a certain threshold calculated from the false alarm probability 
that we can afford we say that the signal is detected. The values of the parameters that maximize A are 
said to be the maximum likelihood (ML) estimators of the parameters of the signal. The magnitude of 
the maximum of A determines the probability of detection of the signal. 

We assume that the noise n in the detector is an additive, stationary, Gaussian, and zero-mean 
continuous random process. Then the data x (if the signal h is present) can be written as 

xit) ^n{t) + h{t). (37) 

The log likelihood function has the form 

logA^{x\h)-^{h\h), (38) 

where the scalar product ( ■ | ■ ) is defined by 

ix\y):=mj^^^^^^df, (39) 

where " denotes the Fourier transform, * is complex conjugation, and Sh is the one-sided spectral density 
of the detector's noise. 

The gravitational-wave signal given by Eq. ( |27| ) consists of two narrowband components around the 
frequencies /o and 2/o and therefore to a very good accuracy the likelihood ratio is given by 

logA = {x\hi) - i(/ii|/ii) + {x\h2) - ^{h2\h2). (40) 



This suggests that we consider the two components of the response function (^7|) as two independent 
signals. Let us take the first component hi of the signal. We can assume that over the bandwidth of the 
signal Shif) is nearly constant and equal to Sh{fa) where fo is the frequency of the signal hi at t — 0. 
Thus in our case the above scalar product can be approximated by 



(x\hi)^^^ I "\it)hiit)dt, (41) 

-To/2 



Shifo) 
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where To is the observation tirae and where the observation interval is [—To/2, To/2]. It is useful to 
introduce the following scalar product 

ix\\y)--=^ [ °' x(i)y(t)dt. (42) 

-To/2 



T 

o 



As long as the detector's noise is stationary over the observation period, this is a good scalar product. 
In realistic observations, the detector's noise will vary slowly during the observation period. We do not 
treat this important issue in this paper. 

The log likelihood function for this signal is approximately given by 



logAi- ^" 



{.x\\h,)-\{h,\\h,) 



(43) 



Sh{h) 

The maximum likelihood estimators can be found by maximizing the following normalized log likelihood 
function 

\ogK[ = (x\\h^)-]^{h,\\h^). (44) 

The normalized log likelihood function does not involve explicitly the spectral density of the noise in the 
detector. 

The signal hi depends linearly on four amplitudes An. The amplitudes depend on the five unknown 
parameters ho, 9, ip, l, and $o and arc independent. The likelihood equations for the amplitudes Axi are 
given by 

^- (-) 

One easily finds that in our case the above set of equations is equivalent to the following set of linear 
algebraic equations 

4 

^M.iiAx^ = {x\\hxi), i=l,...,4, (46) 

where the components M.ij of the 4x4 matrix A1 are given by 

M.,:^{hu\\hii). (47) 

Since over a typical observation time To the phase $ will have very many oscillations, then to a very 
good accuracy we have 

(/iii||M = o, (/iii||/ii4) ^0, (/ii2||M = o, (/ii2||M = o, (48) 

and also 

(/in||/iii) = (/ii3||/ii3) = \A, 

(/112II/112) = (/ii4||/Ji4) = \b, (49) 
(/in 11/112) = (/113II/114) = \C, 
where A := (a||o), B := {b\\b), C := {a\\b). With these approximations the matrix A4 is given by 

^={0 ?)' (^») 
where O is a zero 2 by 2 matrix and C equals 

<^ = Kcb)- pi) 

Thus M. splits into two identical 2x2 matrices. Assuming that a ^ b, A ^ and B ^ the explicit 
expressions for maximum likelihood estimators Au of the amplitudes Au are readily obtained and they 

are given by 

^ ^ BjxWhn) - C{x\\hi2) 
An - 2 , 

2 ^ A{x\\hi2)-C{x\\hn) 

2 ^ Bix\\his)-Cix\\hi,) (52) 
Ais = 2 , 

2 ^ A{x\\hu)-Cix\\hxs) 
Am = 2 , 
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where D is defined by 

D^AB~C^. (53) 
The second partial derivatives of the log likelihood function w.r.t. An are given by 

92 InA'i 



dAudA 



= --My- (54) 



Since a ^ 6 it follows from Schwarz inequality that D > 0. Thus as A > and B > the matrix A4 is 
positive-definite. Therefore the extrema of the log likelihood function w.r.t. An are the local maxima. 
The above ML estimators of the amplitudes An are substituted for the amplitudes An in the likelihood 
function (Q) giving the reduced normalized likelihood function A" = exp(jFi) where !Fi is given by 



BjxWhn? + A{x\\h,2? - 2Cix\\hi,){x\\h,^) 

- D 

BjxWhis)^ + A{x\\hu)^ - 2C{x\\hu)ix\\hu) 
D 



(55) 



Thus to obtain the maximum likelihood estimators of the parameters of the signal one first finds the 
maximum of the functional J^i with respect to the frequency, the spindown parameters and the angles a 
and 5 and then one calculates the estimators of the amplitudes An from the analytic formulae ( |5^ ) with 
the correlations (a;||/iii) evaluated at the values of the parameters obtained by the maximization of the 
functional J^i. Thus we see that filtering for the gravitational- wave signal from a neutron star requires 
four linear filters. Efficient numerical methods to calculate the statistics J^i are discussed in Section 3.4. 

Exactly the same procedure applies to the second component of the signal. The formulae for the 
estimators of the amplitudes A2i and the normalized reduced statistics J^2 are obtained from the above 
formulae by replacing hn by /i2i- 

To consider the optimal detection of the whole two-component signal we need to remember that the 
eight amplitudes An are not independent. They depend on five parameters: ho, 0, ip, t, and <i>o- To find 
the maximum likelihood estimators of the independent five parameters we would have to maximize the 
total likelihood function (given by Eq. (^o|)) with respect to these parameters. This however leads to 
an intractable set of nonlinear algebraic equations which would have to be solved numerically, thereby 
increasing the computational cost of the search for the signal. Instead we propose the following procedure. 

We form the statistics 

^ = S^f' + ^'"^ 

This is just the reduced likelihood function assuming that the eight amplitudes are independent. We 
first maximize the functional with respect to the frequency, spindown parameters and angles a and S 
and we calculate the eight amplitudes from the analytic formulae. We then find the estimators of the 
five independent parameters from the estimators of the amplitudes by least-squares method. We use the 
inverse of Fisher matrix for the covariance matrix in the least-squares method. We shall consider this 
problem in Paper IV. 

To announce the detection of the signal the functional J-' must exceed a certain threshold calculated on 
the basis of the false alarm probability that one can afford. Once T is above the threshold its magnitude 
determines the probability of detection of the signal. Consequently we need to determine the probability 
density function of both when the signal is absent and present. 

We shall first calculate these probabilities when the parameters which !F depends on are known i.e. 
when the filters hu are known functions of time. We shall then explain how to obtain approximate 
formulae for the false alarm and the detection probabilities when parameters of the filters arc unknown. 



3.2 Detection statistics 

We shall first consider the probability density function of the normalized reduced functional J-'i . Let us 

(fc) 

suppose that filters hn are known functions of time, i.e. the phase parameters /o, /o, S are known, 
and let us define the following random variables 

xn ■■= ix\\hn), i = l,...,4. (57) 
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mil 


:= Ei{a;ii} = 


i(AAii 


+ CA12), 


mi2 


Ei{2;i2} = 


i(CAii 


+ Sv4i2), 


mi3 


:= Ei{xi3} = 






mi4 


:= Ei{a:i4} = 




+ Bv4i4). 



Since a; is a Gaussian random process the random variables xu being linear in x are also Gaussian. Let 
Eo{a;ii} and Eijxii} be respectively the means of xu when the signal is absent and when the signal is 
present. One easily gets 

EoixiJ^O, i = l,...,4, (58) 

and 

(59) 
(60) 
(61) 
(62) 

One finds that the covariance matrix for the random variables xu is the same whether the signal is 
present or not and it splits into two identical 2 by 2 covariance matrices C for the pairs (a;ii,a;i2) and 
(xi3, 0:14) of random variables where C is given by Eq. (pT|). Hence the covariance matrix is exactly equal 
to the matrix M given by Eq. ( pO| ) above. 

Thus in effect (xii,xi2) and (0:13, 2:14) are pairs of correlated Gaussian random variables, with pairs 
being independent of each other. For Gaussian variables their first two moments determine uniquely their 
probability density function (pdf). Consequently the joint probability density function p(a;ii, 2:12, Xia, X14) 
is equal to a product Pa{xii, X12) Ph{xiz, XlI^)^ where pa and are bivariate Gaussian pdfs with the same 
covariance matrix C : 

If BS:l^+Axl2-2CS:iiXi2\ 
PQ(a;ii,a;i2) = -^exp I — j, (63) 

and a similar formula for Pbixi3,xi4), where xu = xu when the signal is absent and xu = xu — mu 
when the signal is present. It is interesting to note that when the signal is absent the joint pdf po is 
simply given by 

= -^TF, exp(--?^i): (64) 

where Ti is our optimal statistics. We want to find the pdf of Ti when the signal is absent and present. 
We first decorrelate the variables xu. For the case of Gaussian variables this can always be done by 
means of a linear transformation. Let us consider the following transformation matrix C 

JV o 
o 

where O is a zero 2 by 2 matrix and M is given by 



JV=\ W^±£^ l/^ B + C^ ^^^^ 
-ll^A-C^A/B 1/^B-C^BjA 

Let us introduce new random variables Zu {i — 1,...,4) such that Zu — 'Y^'^j=i^ij^ij- the new 
variables the pdf takes the form 



p(zii, Z12, zi3, Z14) = j—y^ exp 



1 

(2V) 



1 /~2 1-2 1-2 , ~2 \ 
' 2 V^ll "^12 ^13 ' ^14; 



(67) 



Thus zii are independent Gaussian random variables with unit variances. When the signal is absent we 
have Zu = zu and when the signal i 
J^i in the new variables is given by 



have Zu = zu and when the signal is present zu — zu — where = X]t=i ^iji^^ij- The functional 



■^1=0 (^11 + ^12 + ^13 + ^m) ■ (68) 
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The probability density distributions of T\ both when the signal is absent and present are well known. 
When the signal is absent 2T\ has a distribution with 4 degrees of freedom and when signal is present it 
has a noncentral distribution with 4 degrees of freedom and noncentrality parameter A = 

In exactly the same way one obtains the pdf for the normalized reduced functional for the second 
component of the signal. The second component depends on four amplitudes A2i. The decorrelation is 
achieved by the same matrix £. Let us denote the four decorrelated random variables for the second 
component by zn and their means by m!^^ (i = 1, . . . ,4). 

To obtain the pdf of the statistics T for the detection of the two-component signal it is convenient to 
introduce the following normalized random variables z" : 



.H,/^, .r,. = .J^, ^ = l,...,4, (69) 



so that each random variable has a unit variance. Consequently 2T — X^i^il-^D^ has a distribution 
with 8 degrees of freedom when the signal is absent and noncentral x^ with 8 degrees of freedom when 
signal is present. The noncentrality parameter A is given by 

4 4 

After some algebra one finds that A = where 



d ^fW)- (71) 

The quantity d is called the optimal signal-to-noise ratio. It is the maximum signal-to-noise ratio that 
can be achieved for a signal in additive noise with the linear filter ||3l|| . This fact does not depend on the 
statistics of the noise. 

Consequently the pdfs po and pi when respectively the signal is absent and present are given by 

Po{^) = ^exp(-J^), (72) 



pM^) = ^^^^l3{dV2T)exp(^-J^-^d'j, (73) 

where is the modified Bessel function of the first kind and order 3. The false alarm probability Pp is 
the probability that exceeds a certain threshold J-'o when there is no signal. In our case we have 

Pf{^o) ■■= Po(^) dT=[l^To + \tI + exp(-.fo). (74) 

The probability of detection Pd is the probability that T exceeds the threshold To when the signal-to- 
noise ratio is equal to d: 

PD{d,To):= pMT)dT. (75) 

Thus we see that when the noise in the detector is Gaussian and the phase parameters are known the 
probability of detection of the signal depends on a single quantity: the optimal signal-to-noise ratio d. In 
view of its importance we shall investigate in detail the dependence of the optimal signal-to-noise ratio 
on the parameters of the signal in the next section. 

Let us introduce a vector parameter 6^ — (/o, /q, a, 8) denoting the phase parameters. When 
parameters 9^^ are known the optimal statistics T is a, random variable with probability density functions 
given above. When the phase parameters are not known we can think of as a multi-dimensional random 
process J^iO^^) with dimension equal to the number of phase parameters. Such a process is called random 
field. For each realization x{t) of the data random process the corresponding realization of the random field 
is obtained by evaluating J- for filters hu with continuously varying parameters ^f''. For such a process we 
can define the autocorrelation function C just in the same way as we define the autocorrelation function 
for a one parameter random process: 



E[H6f)He\'"')]- (76) 
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Let us first assume that the signal is absent, i.e. x{t) = n{t). For many cases of interest the autocorrelation 
function will tend to zero as the differences = — 9^^ increase. Thus we can divide the parameter 
space into elementary cells such that in each cell C is appreciably different from zero. The realizations 
of the random field within a cell will be correlated (dependent) whereas realizations of the random field 
within each cell and outside the cell are almost uncorrelated (independent). Thus the number of cells 
covering the parameter space estimates the number of independent samples of the random field. For 
some signals the autocorrelation function will depend only on the differences A; and not on the absolute 
values of the parameters. Then the random field T is called a homogeneous random field. In this 
case one can introduce the notion of the correlation hyperellipse as a generalization of the correlation 
time of a stationary process and estimate the area of the elementary cell by the area of the correlation 
hyperellipse. For the general case of a random field the number of elementary cells Nc can be estimated 



from Owen's formula ||29|, 21 with an appropriate choice of the mismatch parameter /i and for the case 
of a homogeneous random field from a formula proposed by one of us p^ . For the parameter values in 
each cell the probability distribution of !F{9^^) can be approximated by probability po(-?^) given by Eq. 
W%- Thus the probability distribution of is given by product of Nc pdfs pq{T). The probability that 



T does not exceed the threshold To in a given cell is 1 — Pf{To), where Pf{To) is given by Eq. (74). The 
probablity that T does not exceed the threshold To in all the Nc cells is [1 — PpiTo)]^" ■ The probability 
Pp that T exceeds To in one or more cell is given by 

P];{To) = l-[l-PF{To)f'. (77) 

This is the false alarm probability when the phase parameters are unknown. When Pf{To) 1 and 
NcPf{To) < 1 we have Pp = NcPf{To)- When the signal is present a precise calculation of the pdf of 
T would be very difficult because the presence of the signal makes the data random process x{t) non- 
stationary. As a first approximation we can approximate the probability of detection of the signal when 
parameters are unknown by the probability of detection when the parameters of the signal are known 
[given by Eq. (^5|)]. This approximation assumes that when the signal is present the true values of the 
phase parameters fall within the cell where T has a maximum. This approximation will be the better 
the higher the signal-to-noise ratio d. An accurate probability of detection can be obtained by numerical 
simulations. Parametric plot of probability of detection vs. probability of false alarm with optimal signal- 
to-noise ratio d as a parameter is called the receiver operating characteristic (ROC). Detailed calculations 
of the number of cells Nc and false alarm probabilities as well as plots of receiver operating characteristic 
for the case of the signal considered here will be given in Paper III. 

The above reasoning is a generalization to the case of many parameters of the idea of an effective 
sampling rate introduced by one of us [p^ and further developed in ||30|] . Related ideas can also be found 
in Ref. ||3^ . 

For large signal-to-noise ratios, the rms errors of the estimators of the parameters of the signal are 
approximately given by the square roots of the diagonal elements of the inverse of the Fisher information 
matrix F with the components given by 



r — 



dh 



(78) 



We shall study these errors in detail in Paper II. For smaller signal-to-noise ratios (e.g. < 10) the errors 
are larger (see Ref. for a discussion in the context of coalescing binaries). 

3.3 Signal-to-noise ratio 

In this subsection we use the following models of the noise spectral densities Sh in the individuals 
detectors. The noise curves for the VIRGO and the initial/ advanced LIGO detectors are taken from 
psf , and the noise curve for the TAMA300 detector is taken from Wideband and narrowband 

versions of the GEO600 detector noise are based on . 

The optimal signal-to- noise ratio d is given by the formula (|7l|): 



d vWO- (79) 

The gravitational-wave signal defined by Eqs. ([l9|)"(^^ consists of two narrowband components around 
the frequencies /o and 2/o and therefore to a very good accuracy the signal-to-noise ratio (|79| ) for that 
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signal can be written as 



dl + dl 



(80) 



where di and d2 are the signal-to-noise ratios for the individual components of the signal. They are given 

by 



di 



To/2 



d2 V{h2\h2) = 

We substitute Eqs. (|l|)-(|2|) to Eqs. (|8^)-(] 
multiples of the frequency /q. We obtain 



[ 2 

I Sfiifo) J-TJ2 



1/2 



[hi{t)Y dt 



To/2 



1/2 



To/2 



[h2(t)fdt\ 



(81) 



(82) 



dl 



dl 



1 

64 



■ sin 2t 



To/2 



To/2 



- (1 -I- cos"' l) 



Fldt 

To/ 2 
-To/2 



and drop out terms which oscillate around zero with 

hi sin^ 261 



1 

16' 



sm^ i 



To/2 

Fldt 

-To/2 
To/ 2 



F^dt -f cos^ i 



Fldt 



To/2 



Shifo) ' 

hi sin* e 

Sh{2fo) ■ 



After performing integrations in (|83|) and (84) we get: 



dl 
dl 



[Ai (5, ^, L) To + Bi (a, <5, V^, T,,)] sin^ 



hi sin^ 26* 



[^2 ((5, ^, l) To + B2 [a, S, L- To)] sin^ 



hi sin* e 



(83) 
(84) 

(85) 
(86) 



The functions Bi and B2 are periodic in the observation time To with the period of two sidereal days (cf. 
Eq. (136) from Appendix B). For simplicity we suppress the explicit dependence of the functions Ak and 
Bk on the angles A and 7. Detailed expressions for the functions and Bk are given in Appendix B. 

For the observation times To longer than several days the signal-to- noise ratios di, d2, and d are 
dominated by terms proportional to the square root of the observation time Tq. This can be seen in 
Figures 2 and 3. 

The signal-to-noise ratios d\ and d\ are complicated functions of the angles a, 5, and 6. We 

have studied the different averages of d\ and d\ over these angles. Averaging is performed according to 
the definition: 



1 

2^ 



2tt 



da X 



1 



dsin(5 X 



1 

2^ 



2tt 



d^p X 



1 



d cos L X 



1 



d0{- 



(87) 



Note that because S e [— 7r/2,7r/2] integration over sin (5 rather than cos^ is involved in Eq. (87). 

Averaging over the angle a discards the oscillatory parts Bi and B2 of the signal-to-noise ratios d^ 



and df : 



{dl)a = Ai (5,V,t)sin\ 



hlTo sin^ 29 
Shifo) ' 



2 ./i^T^sin* 



^^(2/0) 



Further averaging over the orientation angles tp and c gives 



_ 1 , h^IoSU^ 



, 4 . 2 .ft-oToSin^e 

= -e2(^)sm^C " °, . ■ 



(89) 

(90) 
(91) 
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The function 62 (its definition can be found in Appendix B) in the above equations is a fair representation 
of the average sensitivity of a detector at a given location. It depends on the dechnation S of the 
gravitational-wave source as well as on the latitude A of the detector's site and the angle 7 describing the 
orientation of its arms. The product 62 ((5) sin^ C is plotted against the declination 5 for different detectors 
in Figure 4. 

Averaging over the angles a, S, ^, and t yields results which do not depend on the position of the 
detector on the Earth and on the orientation of its arms: 

, ,2^ 1.2. h^To sin^ 29 

Too"-^ 5„(/„) • f''^' 

For the special case of the model of neutron star as a triaxial ellipsoid the angle 9 — 7t/2 and then 
the contribution df to the signal-to-noise ratio vanishes. However for small angles 9 the term may 
dominate over the term ■ The averaging of the above formulae over the angle 9 gives 

We observe that when the noise spectral density at frequencies /o and 2/o is the same the average ( psf ) 
of (i| is more than one order of magnitude greater than the average (|9^) of df . 

We have studied the distribution of the signal-to-noise ratios di, c?2, and d over the angles a, 6, Vi 
and 9 with the aid of the Monte Carlo simulations for the observation time To = 120 days. For each case 
we have generated 10000 sets of angles according to the probability measure defined by the right-hand side 
of Eq. (p7[). We have assumed that the parameter do given by Eq. ( p6| ) is equal to 1. The results are shown 
in Figures 5-7 where we have plotted cumulative distribution functions of the simulated signal-to-noise 
ratios di, ^2, and d for the initial/ advanced Hanford, initial/advanced Livingston, VIRGO, GEO600, and 
TAMA300 detectors. We have performed simulations for two gravitational wave frequencies /q: 100 Hz 
and 500 Hz. The shapes of the distributions of the signal-to- noise ratios di and ^2 do not depend on 
the frequency /o (cf. Eqs. (|8l|)-(^^) and will be the same for nonaxisymmetries generated by different 
physical mechanism e.g. for the case of CFS instability. 

In Table 2 we have given the means and the quartiles for the Monte Carlo simulated cumulative 
distribution functions of the signal-to- noise ratios di, d2, and d for the individual detectors. 

From Figures 5, 6, and Table 2 we see that the simulated distributions of the signal-to-noise ratios 
di, c?2 , and d depend weakly on the position of the detector on the Earth and on the orientation of its 
arms (cf. plots and data for the initial/ advanced Hanford and Livingston detectors). This is related to 



the fact that the averages ( |92| ) and (93) are idependent of the position of the detector on the Earth and 
of the orientation of its arms. 



3.4 Data analysis method 

It is important to calculate the optimum statistics as efficiently as possible. One way to achieve this is 
to take advantage of the speed of the fast Fourier transform (EFT) . Let us consider first the normalized 
reduced functional ^1. One observes that the phase 4> of the signal can be written as (cf. Eq. ([T^)) 

$(t) = 27r/o[t + $„(t; a, S)] + $,(t; a, S), (96) 

where functions $,„ and do not depend on the frequency parameter /q. Let us define the following 
two integrals: 

Fia = I ' 2;(i)a(i)exp[-2$,(i)]exp{-i27r/o[i + $„,(i)]} rfi, (97) 

J-To/2 
rTo/2 

Fib = / 2;(i)6(i)exp[-2$,(i)]exp{-i27r/o[< + $,„(i)]} dt. (98) 

■J-Ta/2 
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detector 




fo = 100 Hz 


/o = 500 Hz 


mean 


90.25 


90.5 




mean 


90.25 


90.5 


90.75 


GEO600 
wideband noise 


di 


0.70 


0.37 


0.72 


1.0 


14. 


7.5 


15. 


21. 




3.1 


0.79 


2.7 


4.7 


21. 


5.4 


18. 


32. 


d 


3.2 


1.1 


2.9 


4.7 


28. 


17. 


28. 


37. 


GEO600 
narrowband noise 


(h 










1.2 


0.()4 


1.3 


1.8 




_ 








180. 


47. 


160. 


280. 


d 


- 


- 


- 




180. 


47. 


160. 


280. 


initial Hanford 


di 


2.8 


1.5 


2.9 


4.2 


48. 


25. 


49. 


70. 


d2 


12. 


3.2 


11. 


19. 


72. 


19. 


63. 


110. 


d 


13. 


4.7 


12. 


19. 


95. 


56. 


96. 


120. 


advanced Hanford 


di 


89. 


46. 


90. 


130. 


480. 


250. 


490. 


700. 




140. 


37. 


120. 


210. 


720. 


190. 


630. 


1100. 


d 


180. 


100. 


180. 


240. 


950. 


560. 


960. 


1200. 


initial Livingston 


di 


2.9 


1.5 


3.0 


4.3 


48. 


25. 


50. 


71. 


d2 


12. 


3.2 


11. 


19. 


72. 


19 


64. 


110. 


d 


13. 


4.8 


12. 


19. 


95. 


58. 


97. 


120. 


advanced Livingston 


di 


89. 


47. 


92. 


130. 


480. 


260. 


500. 


720. 


d2 


140. 


37. 


130. 


220. 


720. 


190. 


640. 


1100. 


d 


180. 


110. 


180. 


240. 


950. 


580. 


970. 


1200. 


VIRGO 


di 


1.5 


0.78 


1.5 


2.2 


46. 


24. 


48. 


68. 


d2 


5.8 


1.5 


5.2 


8.9 


86. 


22. 


76. 


130. 


d 


6.2 


2.3 


5.7 


9.0 


110. 


57. 


100. 


140. 


TAMA300 


di 


0.094 


0.049 


0.098 


0.14 


13. 


6.9 


14. 


19. 


d2 


1.1 


0.28 


0.97 


1.7 


29. 


7.3 


25. 


44. 


d 


1.1 


0.30 


0.98 


1.7 


34. 


17. 


33. 


46. 


initial LIGO/VIRGO 
network 


di 


4.3 


2.3 


4.5 


6.4 


82. 


43. 


86. 


120. 


d2 


19. 


4.8 


17. 


28. 


130. 


35. 


120. 


200. 


d 


20. 


7.1 


18. 


29. 


170. 


100. 


170. 


230. 



Table 2: The means and the quartiles for the Monte Carlo simulated distribution functions of the signal- 

to-noisc ratios rfi , d2, and d for the individual detectors and for the three detector network of the VIRGO 
and two initial LIGO detectors. For the GEO600 detector we use two noise curves: wideband and 
narrowband tuned to 1 kHz with the bandwidth of 30 Hz. We assume that star's ellipticity is 10~^, its 
moment of inertia w.r.t. the rotation axis is 10^'"' g cm^ and its distance from the Earth is 1 kpc. Quantile 
Qx gives a value Zx of random variable z such that probability that z < Zx is less than or equal to x. The 
quantile values at a; = 0.25, 0.5, and 0.75 are called the quartiles. 
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One can write the statistics T\ in terms of the above two integrals as: 

•^1-^ D • (99) 

We can introduce a new time coordinate 

tb{t)=t + <^>„,{t). (100) 

From the expHcit expression for the phase $ given by Eq. ( p^ ) the time shift and its time derivative 
$m can be estimated by 

\^m{t)\ < ^ -5X10^8, 

l*™(t)l < - IX 10-4. 

c 

Assuming the maximum observation time To = 120 days to a very good approximation we have 

n UiTo) = To, -J = 1. (102) 
dtb 

Thus in the new time coordinate the integrals Fia and Fu, can be very well approximated by 

Fia = / °^ x[t(4)]a[t(4)]exp{-i$,[i(ifc)]}exp(-i2^/o4)di6, (103) 

-Fib = f ^ .x[t(tb)]6[t(tb)]exp{-i$,[i(4)]}exp(-i27r/o4)d4. (104) 



Hence we see that with the new time coordinate tb the two integrals ( |103| ) and (|10J) are Fourier trans- 
forms of the functions x[t{tb)]a[t{tb)]exp{—i^s[t{tb)]} and x[t{tb)]b[t{tb)]exp{—i^s[t{tb)]}, respectively. 
To calculate these integrals for a given set of phase parameters we need to perform the following nu- 
merical operations. For the chosen values of the parameters a and S we resample the original time 
series according to the formula ( |100| ) and then wc multiply the resampled time series x{tb) by functions 
a{tb) exp{—i^s[t{tb)]} and 6(tb) exp{— i$s[i(tb)]}. Then we calculate the two Fourier transforms (using 
FFT algorithm). The resampling technique has been proposed by one of us p5[ | and considered as an 
effective data analysis tool for searches of gravitational waves from periodic sources . 
Alternatively one could define new spindown parameters 

A-:-^, fc = l,...,5, (105) 

JO 

and introduce a different time coordinate 

4(i) -t + $„(i) + $;(i), (106) 

where 

j-k+l 



and perform the resampling process according to the formula ( |1G6| ). 

The functions a, b and consequently A, B, C, and D are known and they depend on the declination, 
the right ascension and the time of observation. Their values can be calculated and stored for a fine grid 
of positions of the neutron star on the sky and appropriate observation times before the data analysis is 
carried out. 

The normalized reduced functional for the second component of the signal can be calculated in a 
similar way. Here the corresponding Fourier transforms are given by 

/•To/2 

F2a = / x[titb)]a[titb)]exp{-i2'i>s[titb)]}exp{-iATTfotb)dtb, (108) 

/•To/2 

F2b ^ / x[t{tb)]b[t{tb)]exp{-i2<i>s[t{tb)]}exp{^iAnfotb)dtb. (109) 

J -To/2 
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The statistics T for the whole signal is then calculated from the formula 



T = 



4 



B\F, 



Sh{h)To 



D 



4 



B\F2a? + A|F2b|2 - 2C^[F2aF^^^) 



(110) 



O 



D 



The statistics needs to be calculated on a multidimensional grid of parameter values (excluding the 
frequency parameter /o) covering sufficiently densely the parameter space, and compared against a thresh- 
old. 

4 Networks of detectors 

The analysis of the previous section can be generalized to the case of a network of N interferometers in 
a straightforward manner. Assuming that the noise in each detector is uncorrelated with the others, the 
likelihood function for the network is the sum of the likelihood functions for the individual detectors. 
Therefore we define a statistics J-n for the whole network as the sum of the individual statistics of each 
detector given by Eq. (|5^). We maximize J-n with respect to the phase parameters to obtain their 
estimators. We calculate the estimators of the amplitudes from the analytic fomulae. Then we use a 
least-squares fit to estimate the parameters {ho,9,tp, i-,^o) from the 8A^ amplitude estimators. When 
the phase parameters of the signal are known each of the individual statistics J^i multiplied by a factor 
of 2 has probability density distribution with 8 degrees of freedom when the signal is absent and 
noncental distribution with noncentrality parameter df when the signal is present. The Gaussian 
variables entering each statistics (normalized random variables given by Eqs. (|69|)) have the same 
unit variance. Thus 2^„ has the distribution with SiV degrees of freedom when signal is absent and 
noncentral distribution with noncentrality parameter — when the signal is present. The 

quantity d„ can be defined as the total signal-to-noise ratio of the network. Probability of detection is 
then caclulated by Eq. (^5|). When the phase parameters of the signal are unknown similarly like in the 
one-detector case one can consider a random field which is a sum of the random fields for the individual 
detectors and investigate the correlation function for this random field to obtain independent number of 
cells Nc of the field. One can then calculate the false alarm probability for the network by means of Eq. 



We have studied the distribution of the network signal-to- noise ratios dn2, and the total network 
signal-to-noise ratio dn = \/ rf^i + c?^2 '^^'^^ ^^le angles a, S, t/j, l, and 9 with the aid of the Monte Carlo 
simulations for the observation time To = 120 days. We have restricted ourselves to the three detector 
network of the VIRGO and two initial LIGO detectors. For each case we have generated 10000 sets of 
angles according to the probability measure defined by the right-hand side of Eq. (|^ . We have assumed 
that the parameter do given by Eq. (^) is equal to 1. The results are shown in Figure 8. We have 
performed simulations for two gravitational wave frequencies /q: 100 Hz and 500 Hz. 

In Table 2 we have put the means and the quartiles for the Monte Carlo simulated cumulative 
distribution functions of the signal-to-noise ratios dni, (in2, and c?„ for the three detector network of 
the VIRGO and two initial LIGO detectors. Adding the GEO600 and TAMA300 detectors wiU not 
significantly change these signal-to- noise ratio values, but the smaller detectors can play an important 
role in making coincident detections by improving the confidence that the candidate events registered by 
larger detectors are not due to un-modelled noise. 
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A The phase of the gravitational- wave signal 



We assume that in the rest frame of the neutron star the time dependence of the phase of the gravitational- 
wave signal can be written as a power series of the form: 



rk+l 



k=0 



(fc + 1)!' 



(Ill) 



where r is the proper time in the neutron star rest frame. The assumption (111) means that the instan- 
taneous frequency of the signal in the rest frame of the neutron star is given as 



/„s(t) 



1 rf^ns(r) 
2tt dr 



-V - 



fc=0 



(112) 



(k) 

so /ns is the kth time derivative of the frequency evaluated at t = 0. 

We assume that the neutron star is moving with respect to the SSB uniformly along a straight line 
according to the equation 

rns{t) ^ roUo + VnsUy (^t + , (113) 

where ro ;= |rns(i = ^^o/c)|, no :— r^sit ~ — ro/c)/ro. If we denote by the constant velocity vector 
of the neutron star then v^s ■= |vns| and rii, := Vns/vns- The time t in Eq. (113) is the time coordinate 
in the SSB rest frame. We do not allow the neutron star to have an intrinsic acceleration. This means 
we exclude binary neutron stars, except of the binary periods so long that the acceleration effects may 
be accurately approximated by a Taylor series during the observation time. 

The phase observed at the SSB at some time t was emitted by the star at the coordinate time t' such 
that 

\rns{t')\ 



c 

One can show that the relation between the time t' and the star's proper time t is as follows 

rQ_ 

c 



(114) 



(115) 



where /3ns ■— Vns/c. In Eq. (115) the time dilation effect is taken into account. We have also assumed 
that T — when the star's position vector w.r.t. the SSB is rong. We can now write 



SSB 



{t) = *ns(r(i)), 



(116) 



where 5'ssb(^) is the phase observed at the SSB at time t, and the time t can be expressed in terms of t 
by means of Eqs. (114) and (115). 



Collecting Eqs. (Ill) and (114)-(11(;) together we can write 



*SSB(t) ^^0 + TT^ (1 - Pnsf^'^^" (t'it) 



c 



fe+i 



where t' is the solution of Eq. (114) for a given time t. It reads 
t' 



t -f —(3ns [/3ns + (no • n„) 



I3U^ + 2^/3„, [A.S + (no • n,)] < + ^ [1 + /5ns (no • n,) 



(117) 



(118) 



We expand the function ^'ssb given by Eqs. ( [117|) and ( |118| ) w.r.t. time t around t = Q. The first few 
terms of the expansion read 



^'ssB(i)-$o _ <°> , 

" — /sSB * 



27r 



/sSB 



(0) Y_ 

(l + (no-n„)/3ns)2(ro/c) ^^^^] 2 



\ (2) 

+\.h 



SSB 



3 ((no ■ n.)^ - 1) (31, 
(1 + (no • n^)/3ns)^ (r-o/c) 



(1) 

/sSB — 



(/3ns + (no ■ n^))/3ns ' 

(1 -I- (no • n„)/3ns)^ (ro/c) 



-+0(t^), (119) 
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where 



(fe) 

/sSB 



,2 \ 



(1 + (no • n„) /3ns) 



fe+ 



]^ /ns; k — 0, . . . , ; 



(120) 



As a result of the motion of the neutron star w.r.t. the SSB th e Tay lor expansion ( |119[) of the phase 
contains infinitely many terms, even if we restrict, as in Eq. (|l 1 1| ) , the intrinsic spindown of the star to 
finite number of terms. When the neutron star moves radially w.r.t. the SSB then (no • n„) ~ 1 and the 
function 'I'ssB can exactly be written as the finite sum: 



*ssB(t) = $0 + 27r }^ /gg3___ 



(121) 



fc=0 



We shall assume the following polynomial model of the phase of the gravitational radiation observed 
at the SSB: 

*SSB(i) = 'i>0 + 2^^ /oy^— ry, (122) 



where the new spindown parameters /q do not in general coincide with the Doppler scaled intrinsic 



spindown parameters /sge defined by Eq. (|120D . 

We write the position vector rd of the detector with respect to the SSB as 

rd(i) = rd(i)nd(i). 



(123) 



where r^it) '■— |rd(OI ^-^^d nd(t) Yd{t)/ri{t). The phase of the gravitational- wave signal at the time t 
at the detector's location corresponds to the phase near the neutron star at an earlier instant of time <", 
where t" is the solution of the equation 



t = t 



n , |rns(t")-rd(i)| 



(124) 



The same value of the phase is observed at the SSB at time 

^, , |rn.(t-)| 



thus using Eq. ( |122| ) we can write 



fe=0 



(fc 



|r„s(t"(0)l 



fc+i 



(125) 



where t" (t) is the solution of Eq. (124) for given time t. Using Eqs. (113) and ( |123| ) we express the 
solution t" to Eq. (124) in terms of the time t and the two small parameters /3ns and x :— r^jri^: 



t"{x,(3ns) = Y 



0'^ 



t H /3ns [(no • n„) - (nd • n„)a; + /3ns] 



t + ^/3ns [(no • n„) - (nd • n^)x + /3ns] j + (l - PL) 
-| [1 + I3l^ - 2(no • nd)a; + x^ + 2 ((no • n„) - (nd • n„)a;) /3„s] - 



1/2 



(126) 



Using Eq. (123) we also find that 



|rns(t")l 



? + + (no-n.)/3ns+(t" + ^) A^. 



(127) 



We now study how to simplify the phase '^^i given by Eqs. (125)-(127). 
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An optimal method to detect our signal in noise developed in Section 3 involves correlating the data 
with templates of the signal. In general if the phase of the template differs from that of the signal 
by as little as 1/4 of a cycle the correlation will be significantly reduced. Thus we adopt the criterion 
that we exclude an effect from the model of the signal in the case when it contributes less than 1/4 of 
a cycle to the phase of the signal during the observation time. That this criterion is only a sufficient 
condition but not necessary follows from the correlations among parameters of the phase. The shifts in 
the values of the parameters in the template phase away from the true values of the parameters in the 
signal phase can compensate for the effects in the signal not taken into account in the templates. This 
effect was observed for the case of coalescing binaries |3|, ^ . Finally we stress that such shifts 

in the template parameter values mean that the estimators of the parameters of the signal when using 
an inaccurate template will be biased. It may happen that these biases are much larger than the rms 
errors of the estimators. Thus templates accurate to 1/4 of a cycle over the observation time may not 
be needed to detect the signal, but they will be needed to obtain accurate estimates of the errors in 
parameter measurements. 

In calculating the number of cycles we assume a long observation time of 120 days, the maximum 
gravitational wave frequency of 1 kHz, and the extreme case of a neutron star at a distance rg = 40 pc 
with Vns = 10"^ km/s. For this extreme case the parameters x and /3„s assume the values (as to a good 
approximation rd = 1 AU): 

X = 1.21 X 10"^ /3„s = 3.34 X 10"^ (128) 
The numerical values of the spindown parameters /q we estimate by means of the relation: 



fciA, (129) 



where /o is the radiation frequency and r is the spindown age of the neutron star. As the extreme case 
we will consider r = 40 years. 

It is convenient to carry out the Taylor expansion of the phase ( p_25[ ) with respect to the parameters 
X and /3ns- We note that for any n 

^(. = 0,/.„. = 0)^0. 

Analysis of the first few terms of the Taylor expansion shows that for the observation times Tq < 120 
days, neutron star distances tq > 40 pc, velocities Vns ^ 10^ km/s, frequencies /o < 1 kHz, and spindown 
ages r > 40 years, the only terms which can contribute more than 1/4 of a cycle to the phase of the 
signal, read 

4^ (fc) fk+l On f V , \ (fc) 

fc=0 ^ ' > \ u / 

where Vnsj^ := Vns— (no • Vns) no is the component of the velocity Vns perpendicular to the vector no. The 
ratio Viisi/^'o determines the proper motion of the star. The term in the above expansion proportional 
to v„s_l/7'o contributes at most ~4 cycles. We shall not consider it in this paper. We shall look at the 
possibility of its determination in the next paper of this series. Consequently we restrict ourselves to a 
phase model at the detector of the form 

^ ^ (fc) fk+x ^ (fc) ^fe 

vfd = $0 + /o7^-W + -no • rd ^ k^v (131) 

fc=0 ^ '' k=0 



The model (131) contains the position rd of the Earth relative to the SSB, which we now consider. In 



addition we must consider extra, purely relativistic effects left out of (131). 

Motion of the Earth w.r.t. the SSB is very well determined and there are several computer ephemeris 
routines available In this paper we assume for simplicity that the Earth moves on a circular orbit 
around the Sun. The eccentricity of the Earth's orbit (e® = 0.017) introduces a change of about 8.3 x 10"^ 
cycles in the phase w.r.t. the phase for circular orbit for 1 kHz signal, so it must be included in realistic 
filters. But it introduces no new parameters so we ignore it here. We also ignore the motion of the Earth 
around the Earth-Moon barycenter. 
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There are two types of relativistic corrections. One originates in the difference between the coordinate 
time t which we used in the derivation of the phase model and the proper time r in the detector's reference 
frame. The difference is due to the combined effect of the gravitational redshift and the time dilation. 
The other correction is the Shapiro delay caused by propagation of the gravitational wave through the 
curved spacetime of the solar system. We estimate the contribution to the number of cylces in the phase 
produced by these corrections. 

The difference between the coordinate time t in the first order post-Newtonian coordinate system 
which is assumed to be the rest frame of the SSB and the proper time r kept by a terrestial clock is 
discussed in detail in Ref. Ell]. The difference := t — r is given by the integral 



1 



72 iU[rit')] 



<t') 



dt' 



(132) 



where r is the position vector of the clock w.r.t. the SSB, v := r is the clock's coordinate velocity, and 
U [r(t)] is the instantaneous gravitational potential at the clock's location. The time difference described 
by the integral (132) can be split into the secular and periodic part. The secular difference is due to the 
practically constant rotational velocity and the Earth's gravitational potential at the detector's location 
as well as the average orbital velocity of the Earth and the average gravitational potential along the 
Earth's orbit. This secular difference corresponds to the rescaling of the time coordinate and can be 
incorporated in to th e definition of the spindown parameters. The main contribution to the periodic part 
of the integral (|l32| ) was calculated by Clemence and Szebehely and then corrected by Blandford and 
Teukolsky [H. It can be written as 



2GMfT,e, 



periodic 



(1-4) 



1 



1 



sin M® ■ 



1 . 3 2 

-em sm2MfB + -6^ 



sin iMa 



(133) 



where Mq is the mass of the Sun, a® — 1 AU, fJo is the mean orbital angular velocity of the Earth, 
and M® are the eccentricity and mean anomaly of the Earth's orbit. The quantity (A^;" 



periodic 



varies m 



time with the period of one year and has the amplitude ~ 1.7 x 10"'^ s, so for a 1 kHz gravitational wave 
the contribution of this correction to the total number of cycles is not greater than ~2 cycles. Even when 
it must be included in a filter, it introduces no new parameters. 

The magnitude of the Shapiro delay can be estimated from the relation jii) (neglecting the eccentricity 
of the Earth's orbit) 



A. 



2GMc. 



log 



1 



1 



cos ( 



(134) 



where 9 is the star-Sun-detector angle at the time of observation. To estimate the maximum value of the 
Shapiro delay we consider a neutron star in such position that at some instant of time the line of sight 
from the detector to the neutron star is tangent to the surface of the Sun. Then 6* = 6*1 ~ tt — C, where 
Q ~ Rq/1 AU ~ 4.65 X 10"'^ rad {Rq is the radius of the Sun). Six months later 6* = 6*2 — Cj so the 
amplitude of the correction is 



2GMq , 1 
5— log 



cos 02 



1 + COS 9i 



1.2 X 10" 



For a 1 kHz gravitational wave this gives ~0.1 cycles. So the Shapiro delay will be unobservable. 

We see that the relativistic corrections that need to be applied to our formula are small. By our 1 /4 
of a cycle criterion they can be neglected if we search for signals with frequencies less than '^lOO Hz. We 
shall not consider these corrections in this and the following papers of the series since they are unlikely 
influence our results. However they may need to be included in filters. 

B Signal-to-noise ratio 

In this appendix we give the detailed expressions for the functions Ai, A2, and B2 from Eqs. ( p5| ) 
and (H). They read: 



Ak (5, ip, i) 



Fk{L)ei{5) coslV- + Gfe(0e2((5), 



(135) 
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1 ^ / \ 

Bk{a,5,ip,L]To) = ^ X! ( ) 

^ n— 1 ^ 



X {Ckn (<5, V', COS [n{a - 0^)] + -Dfc„ {5, ip: l) sin [n(a - (t>r)]} ■ 
where fir is the rotational frequency of the Earth, To is the observation time, and where 
Ckn{S,ip,i.) = Fk{L)[fi„{6)cos4:tl; + gin{S)sm4:tl;]+Gk{i.)hin{S), 

Dkn{5,'4),L) = Fk{i)[f2n{5) cos Alp + g2n{6) sin -il};]^ Gk{i)h2n{&)-, 



Fi{i) = -— sin*t. 



F2{i) 



■ 4 

■ sm t, 



Gi(t) = ^sin^t(l + cos^t) , G2(t) = ^ (l + 6cos2 i + cos'' t) 



(136) 

(137) 
(138) 

(139) 
(140) 



The functions ei, 62, fkn, and (7fc„ (fc = 1, 2, n = 1, . . . , 4) entering Eqs. (135), (137), and ( |138| ) are equal 
to 



ei(5) 


= 4ji cos^ 5, 


62(5) 


= 4j2 - is cos 2(5 + ji cos^ 2(5, 




= —4^4 cos'^ 5 sin (5, 


/l2('5) 


= cos^ 5 (3 — cos 2(5) , 


/l3(<5) 


= — je (7 — cos 2(5) sin 25, 


/l4('5) 


= -jy (35 - 28 cos 2(5 + cos 45), 


/2l(<5) 


= — 2 8 js cos^ (5 sin 5, 


/22(<5) 


= —7jr) (3 — COS 25) cos^ 5, 


/23(^) 


= ~ jio (7 — cos 25) sin 25, 


/24('5) 


= -ill (35- 28 cos 25 + cos 45) , 
= 28^9 cos^ 5 sin 5, 


511 (5) 


= 28j8 cos^ 5, 


512(5) 


513 (<5) 


= 2jio (5 — 3 cos 25) cos 5, 


514(5) 


= 16jii (3 — cos 25) sin 5, 


521(5) 


= —4^4 cos'^ 5, 


522(5) 


= 4j5 cos^ 5 sin 5, 


523(5) 


= — 2 je (5 — 3 cos 2(5) cos (5, 


524(5) 


= — IGjV (3 — cos25) sin5. 


/lll(5) 


= (ji2 - j4C0s2(5)sin25, 


hi2i5) 


= (ji3 - 35 cos 25) cos^ 5, 




= 4^6 cos^ 5 sin (5, 


hi4{5) 


= —8^7 cos'' 5, 


h2i{5) 


= js (1 - 7 cos 2(5) sin 25, 


^22(5) 


= "Jg (5—7 cos 25) cos^ 5, 


h2M 


= AjiQ cos^ 5 sin 5, 


^24(5) 


= ~8iii cos^ 5, 



where the coefficients ji , 



ji (A, 7) 
i2 (A, 7) 
js (A, 7) 
J4 (A, 7) 
is (A, 7) 
je (A, 7) 
J7 (A, 7) 
is (A, 7) 
ig (A, 7) 

110 (A, 7) 

111 (A, 7) 

112 (A, 7) 

113 (A, 7) 



ji3 depend on the angles A and 7: 
1 



256 
1 



(4 - 20 cos^ A + 35 sin^ 27 cos'' A) , 

- (68 - 20 cos^ A - 13 sin^ 27 cos" A) , 
1024 ^ in 

(28 - 44 cos^ A + 5 sin^ 27 cos"* A) , 

128 

— (2-7sin^ 27Cos2 A) sin 2A, 

^ (3 - 7 COS 47 - 7 sin^ 27 cos^ A) cos^ A, 

— (2 cos 47 + sin^ 27 cos^ A) sin 2A, 
96 ^ ' 



1 



(4 cos 47 sin^ A — sin^ 27 cos" A) 



1024 

1 „ 

— sin47COS A, 

1 „ 

— sin 47 cos A sin A, 
32 ' 

—— sin 47 (5 — 3 cos 2A) cos A, 
192 

sin 47 (3 — cos 2A) sin A, 

1024 ' ' 

— (l4-sin^ 27Cos^ A) sin 2A, 
32 

7^ (9 — 5 cos 47 — 5 siii^ 27 cos^ A) cos^ A. 
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C Monte Carlo analysis of the two signal-to-noise ratio compo- 
nents 

We have also studied the distribution of the individual signal-to-noise ratios and d"^ given by Eqs. ( p5| ) 
and ( p6| ) over the angles a, (5, ip, and t by means of the Monte Carlo simulations (for the observation 
time To — 120 days). The results are shown in Figures 9-13. The signal-to-noise ratios d^ and ci| are 
normalized here by means of the quantities 

('^l(120days)) (^1 ) a.5,iA,<. |T„ = 120days , (c^2(120days)) •= {d-l) a,S,ii,L\Ta = 120days- 

In Tables 3 and 4 we summarize the statistical characteristics of the simulated distributions of the 
normalized signal-to- noise ratios df and d"^- We give extremal values, means, standard deviations (std), 
and medians. 



detector 


min 


max 


mean 


std 


median 


GEO600 


0.0 


2.0 


1.0 


0.45 


0.96 


LIGO Hanford 


0.0 


1.7 


1.0 


0.41 


1.0 


LIGO Livingston 


0.0 


1.5 


1.0 


0.34 


1.1 


VIRGO 


0.0 


1.6 


1.0 


0.36 


1.1 


TAMA300 


0.0 


1.9 


1.0 


0.38 


1.0 



Table 3: Statistics of d?/(d?(i20days))- 



detector 


min 


max 


mean 


std 


median 


GEO600 


0.18 


4.0 


1.0 


0.72 


0.77 


LIGO Hanford 


0.12 


3.3 


1.0 


0.68 


0.79 


LIGO Livingston 


0.27 


2.7 


1.0 


0.64 


0.80 


VIRGO 


0.26 


3.2 


1.0 


0.66 


0.79 


TAMA300 


0.18 


2.8 


1.0 


0.64 


0.80 



Table 4: Statistics of c'2/('^2(i20days))- 
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Figure 1: Power spectrum of the noise- free response of an interferometer located near Hannover to 
gravitational-wave signal from the Crab pulsar at twice the rotation frequency. We have assumed the 

(1) (2) 

frequency /o = 29.937 Hz and the spindown parameters /q = -3.773 x 10"^" s"^ /„ = 0.976 x 10"^° s~^ 

(3) 

/q = —0.615 X 10~^° s~^. A 24-day long signal was analysed corresponding to the frequency resolution 
of around 4.8 x 10~^ Hz. The power spectrum shows 24 main peaks resulting from the periodic phase 
modulation of the signal. In each interval between the main peaks there are additional subsidiary peaks 
arising from the amplitude modulation of the signal. 
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Figure 2: The relative contribution of the part dgqr of the signal-to-noise ratio proportional to the square 
root of the observation time To to the total signal-to-noise ratio d for the GEO600 detector with the 
wideband noise curve (4qr ~ hoSinC[AisiT? 20 / Sh{fo) + A2iiin^ 6 / Sh[2fo)Y/'^^/To, cf. Eqs. (||), 
and (p6|)). A hypothetical neutron star is assumed to be in the distance of 40 pc from the Earth and 
to emit gravitational waves with frequency /o = 100 Hz, star's ellipticity is 10^^, its moment of inertia 
w.r.t. the rotation axis is lO'*^ g cm^. We also set a - 0^ = 15°, 5 = 35°, ip = 11.25°, l = 22.5°, and 
9 = 45°. 
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Figure 3: The total signal-to- noise ratio d as a function of the observation time To for the GEO600 
detector (with the wideband noise curve). The neutron star parameters are the same as in Figure 2. 
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Figure 4: The plot of the product e2{5) sir? against the declination 6 of the gravitational- wave source 
for different detectors. It can be shown (of. Appendix B) that for a hypothetical detector located at the 
latitude A = ± arccos a/2/3 w ±35.26° there exist eight different orientations 7 of its arms such that the 
function 62 = 1/5, i.e. 62 does not depend on the declination 5 of the gravitational-wave source. These 
orientation angles are: 70, 90° ± 70, 180° ± 70, 270° ± 70, 360° - 70, where 70 = | arcsin -y/s/S w 25.38°. 
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Figure 5: Cumulative distribution functions of the simulated signal-to-noise ratios di, d2, and d for the 
VIRGO (V), initial Hanford (H), and initial Livingston (L) detectors. Wc assume that star's cUipticity 
is 10~^, its moment of inertia w.r.t. the rotation axis is 10^^ g cm^ and its distance from the Earth is 1 
kpc. The observation time is 120 days. The left column is for /o = 100 Hz and the right one is for /o = 
500 Hz. Wc have also shown the cumulative distribution function of the signal-to- noise ratio d2 for the 
GEO600 (G) detector with the narrowband noise tuned to 1 kHz with the bandwidth of 30 Hz. 
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Figure 6: Cumulative distribution functions of the simulated signal-to- noise ratios di, d2, and d for the 
advanced Hanford (H) and advanced Livingston (L) detectors. The observation time is 120 days. The 
left column is for /o = 100 Hz and the right one is for /o = 500 Hz. The neutron star parameters are the 
same as in Figure 5. 
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Figure 7: Cumulative distribution functions of the simulated signal-to-noise ratios di, d2, and d for the 
GEO600 (G) with the wideband noise curve and TAMA300 (T) detectors. The observation time is 120 
days. The left column is for /o = 100 Hz and the right one is for /o = 500 Hz. The neutron star 
parameters are the same as in Figure 5. 
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Figure 8: Cumulative distribution functions of the simulated signal-to-noise ratios rf„i, d„2, and rf„ for the 
three detector network of the VIRGO and two initial LIGO detectors. The observation time is 120 days. 
The left column is for /o = 100 Hz and the right one is for /o = 500 Hz. The neutron star parameters 
are the same as in Figure 5. 
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Figure 9: Histograms of the simulated probability density and cumulative distribution functions of the 
normalized signal-to-noise ratios c^i/(c^i(i20days)) ^^'^ '^2/('^2(i20days)) GEO600 detector. 
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Figure 10: Histograms of the simulated probability density and cumulative distribution functions of the 
normalized signal-to-noise ratios c^i/(c^i(i20days)) ^^'^ '^2/('^2(i20days)) LIGO Hanford detector. 
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Figure 11: Histograms of the simulated probability density and cumulative distribution functions of the 
normalized signal-to-noise ratios c^i/(c^i(i20days)) ^^'^ '^2/('^2(i20days)) LIGO Livingston detector. 
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Figure 12: Histograms of the simulated probability density and cumulative distribution functions of the 
normalized signal-to-noise ratios c^i/(c^i(i20days)) ^^'^ '^2/('^2(i20days)) VIRGO detector. 
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Figure 13: Histograms of the simulated probability density and cumulative distribution functions of the 
normalized signal-to-noise ratios c^i/(c^i(i20days)) ^^'^ '^2/('^2(i20days)) TAMA300 detector. 
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